Search for correlation between GRB's detected by BeppoSAX 
and gravitational wave detectors EXPLORER and NAUTILUS 
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Abstract 

Data obtained during five months of 2001 with the gravitational wave (GW) detec- 
tors EXPLORER and NAUTILUS were studied in correlation with the gamma ray 
burst data (GRB) obtained with the BeppoSAX satellite. During this period Bep- 
poSAX was the only GRB satellite in operation, while EXPLORER and NAUTILUS 
were the only GW detectors in operation. No correlation between the GW data and 
the GRB bursts was found. The analysis, performed over 47 GRB's, excludes the 
presence of signals of amplitude h > 1.2 X 10~ 18 , with 95% probability, if we allow a 
time delay between GW bursts and GRB within ±400 s, and h > 6.5 x 10 -19 , if the 
time delay is within ±5 s. The result is also provided in form of scaled likelihood 
for unbiased interpretation and easier use for further analysis. 
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1 Introduction 

One of the most important astrophysical phenomena still lacking an explanation is the 
origin of the celestial gamma-ray bursts (GRB). These are powerful flashes of gamma- 
rays lasting from less than one second to tens of seconds, with isotropic distribution in the 
sky. They are observed above the terrestrial atmosphere with X- gamma-ray detectors 
aboard satellites |], ||]. Thanks to the BeppoSAX satellite 0, afterglow emission at lower 
wavelengths has been discovered f|, |j| U and we now know that at least long (> Is) 
GRB's are at cosmological distances, with measured red shifts up to 4.5 (see, e.g., review 
by Djorgovski J7| and references therein). Among the possible explanations of these events, 
which involve huge energy releases (up to 10 54 erg, assuming isotropic emission), the most 
likely candidates are the collapse of a very massive star (hypernova) and the coalescence of 
one compact binary system (see, e.g., reviews by Piran M and Meszaros || and references 
therein). In both cases emission of gravitational waves (GW) is expected to be associated 
with them (e.g. Ref. |I(J). According to several models, the duration of a GW burst is 
predicted to be of the order of a few milliseconds for a variety of sources, including the 
coalescing and merging black holes and/or neutron star binaries. Therefore GW bursts 
can be detected by the present resonant detectors, designed to detect GW through the 
excitation of the quadrupole modes of massive cylinders, resonating at frequencies near 1 
kHz. 

At the distances of the GRB sources (~ 1 Gpc), the GW burst associated with a 
total conversion of 1-2 solar masses should have amplitude of the order of h ~ 3 x 10~ 22 . 
The present sensitivity for 1 ms GW pulses of the best GW antennas with signal to 
noise ratio (SNR) equal to unity is h ~ 4 x 10~ 19 (see e.g. Ref. [|TTJ), which requires a 
total conversion of one million solar masses at 1 Gpc. However, although detection of a 
gravitational signal associated with a single GRB appears hopeless, detection of a signal 
associated with the sum of many events could be more realistic. Thus we launched a 
program devoted to studying the presence of correlations between GRB events detected 
with BeppoSAX and the output signals from gravitational antennas NAUTILUS and 
EXPLORER. 

Searching for correlation between GRB and GW signals means dealing with the 
difference between the emission times for the two types of phenomena. Furthermore, 
there is also the fact to consider that the time difference can vary from burst to burst. 
In the present analysis we use an algorithm based on cross-correlating the outputs of two 
GW detectors (see JL2], |i~3"[), thus coping with the problem of the unknown possible time 
difference between GRB and GW bursts, and also of the unmodelled noise. 



Experimental data 



The Rome group operates two resonant bar detectors: EXPLORER [14]], since 1990, at the 
CERN laboratories, and NAUTILUS since 1995, at the INFN laboratories in Frascati. 
The two detectors, oriented nearly parallel, are very similar. They have operated over the 
last few years with various levels of sensitivity, and since March 1 st 2001 they have been in 
operation simultaneously with the best ever reached sensitivity for millisecond bursts, of 
the order of h « (4 — 5) x 10~ 19 . The detectors consist of massive cylindrical bars 3 m long 
made of high quality factor Aluminum alloy 5056. The GW excites the first longitudinal 
mode of the bar, which is cooled to liquid helium temperature to reduce the thermal noise. 
To measure the bar strain induced by a GW, a secondary mechanical oscillator tuned to 
the cited mode is mounted on one bar face and a sensor measures the displacement of the 
secondary oscillator. 
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Table 1: Main characteristics of the two detectors. / indicates, for each detector, the two 
resonant frequencies and Af indicates the bandwidth. The relatively larger bandwidth of 
EXPLORER is due to an improved readout system. 



detector 


latitude 


longitude 


orientation 


mass 


/ 


T 


A/ 










kg 


Hz 


K 


Hz 


EXPLORER 


46.45° N 


6.20° E 


39° E 


2270 


904.7 


2.6 


« 9 












921.3 






NAUTILUS 


41.82° N 


12.67° E 


44° E 


2270 


906.97 


1.5 


» 0.4 












922.46 







The data have a sampling time of 4.544 ms and are processed with a filter matched 
to delta-like signals for the detection of short bursts |LlJ . The filter is adaptive and 
makes use of power spectra obtained during periods of two hours. The filtered output 
is squared and normalized using the detector calibration, such that its square gives the 
energy innovation x(t) for each sample. In the presence only of well behaved noise due to 
the thermal motion of the bar and to the electronic noise of the amplifier, the probability 
density function of x(t) is 



f(x) = = exp 

V 27rT eff X 



X 



2T, 



eff 



where x(t) is expressed in kelvin units, and the average value of x, T e ff, called effec- 
tive temperature, gives an estimation of the noise. If a signal of energy E due to an 
impulsive force acting on the bar is generated at time t a , the change of the filtered data 
energy with time, neglecting the noise contribution, has an envelope which depends on 
the detector bandwidth as follows: 

E s {t) — E ■ exp [— 2-7T \t — 1 \ Af] (2) 

where the bandwidth Af is given in Table [I]. In addition to the well behaved and mod- 
eled noise (electronic and thermal noise), other sources of noise are active, sometimes of 
unknown origin. This requires more than one detector to be simultaneously used, in order 
to discriminate a real signal from noise. 

For our detectors, the relationship between burst energy E expressed in kelvin and 
dimensionless amplitude h is given by [[L7j 



E=7^^^A^-) 2 m (3) 



(7.97 x 10- 18 ) 2 VI ms 

where r gw is the duration of the burst, conventionally assumed to be r gw = 0.001 s. 

In the following, for a given signal of energy E, we shall use the signal-to-noise ratio 

SNR = JL. (4) 

For the GRB's we consider the observations made with the Gamma Ray Burst 
Monitor (GRBM, fL9|) aboard the BeppoSAX satellite, the only satellite for the GRB 
detection in operation during 2001. The GRBM is an all sky monitor which operates in 
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Figure 1: Distribution of T e ff [K], in the period 1 March through 17 July 2001. 



the range from 40 to 700 keV, with a GRB detection rate of about 0.7 events/day. For 
each GRB which triggers the on-board trigger logic, high time resolution (up to 0.5 ms) 
time profiles are transmitted. In addition, 1 s count rates from the GRBM in two energy 
bands (40 — 700 keV and > 100 keV) are continuously recorded and transmitted. From 
the GRBM data rough information on the GRB direction can also be derived p0| . 

For the present analysis we use two quantities from the GRBM data: the initial time 
of each GRB and its burst duration. In the period 1 March through 17 July 2001 we have 
101 BeppoSAX bursts, but only 51 occur at times when both antennas were operating. 
The GRB fall into two categories: one of the 38 GRB events with short duration (< Is) 
and another one of 63 GRB with longer duration. The short duration events include all 
the GRB's which did not trigger the on-board logic, thus their exact duration is not 
available. Indicating with t 7 (trigger time) the initial time for the GRB's, for each burst 
we considered the EXPLORER and NAUTILUS data in t m = 800 s intervals centered 
at the try (t,y ± 400 s). Each interval is covered by 800/0.004544 = 176056 data samples. 
The average value of the data gives, in absence of signal, noise temperature T e ff. The 
distributions of T e ff obtained for these stretches of GW data are shown in Fig. |l|. We 
note the good behavior of both detectors. In view of this good behaviour we thought to 
limit our analysis to the 47 stretches with T e ff < 20 mK on both detectors, losing only 
four stretches. For these stretches we have the following average values: T e e X ff = 8.8 mK 
and T%ff = 6.1 mK. 



3 Cross-correlation analysis 

Due to the unknown time gap between GW's and GRB's, it is essential to have at least 
two independent GW detectors at disposal. In such a case we can cross-correlate the 
two outputs, since we expect that, if any GW burst arrives on the Earth, both detectors 
respond at the same time within the travel time difference of the GW burst. EXPLORER 
and NAUTILUS being about 700 km apart, the maximum travel time difference is about 
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Figure 2: The cross-correlation R(t), averaged over the 47 GRB's versus the time shift r 
in seconds. No positive correlation at r = s is visible. 

3 ms, shorter than the time uncertainty of the measurements, of the order of 1 sampling 
time. Since we use resonant transducers, the signal is, for each detector, distributed over 
the two resonance modes, with a beat period of 64 ms. To cope with this problem we 
averaged our data over 16 samples, that is over a time of 16 x 0.004544 = 0.0727 s. Thus 
during each 800 s interval we have 11004 pairs of data which we can cross-correlate. 
The correlation function is defined by 



where the summations are extended over the 11004 pairs, x(U) refers to the EXPLORER 
squared data and y(ti) is the same quantity for NAUTILUS. r(r) is dimensionless, by 
definition. If simultaneous signals due to GW bursts arise both in the EXPLORER and 
NAUTILUS detectors, no matter when, with respect to the GRB arrival time, but within 
the time window of ±400 s, we should find a larger value of r(r) for r = s. 

The 47 values of r(r), one for each GRB, are used to calculate the average cross- 
correlation R(t), shown in Fig. |2|. No significant correlation is visible. In particular, at 
r = 0s we get the negative value of —0.0027. We also repeated the above calculations 
separately for the GRB's with short (< 1 s) and long (> 1 s) time duration. Again, no 
correlation appears between EXPLORER and NAUTILUS using a window of ±400 s 
centered at the GRB arrival times. 

3.1 Modeling the average cross-correlation 

In order to provide quantitative information out of this null result, we need to model our 
expectations for the average cross-correlation at zero delay time R(0), hereafter indicated 
with R , under the hypothesis that such signals do exist. 

This modeling means going through the following considerations. 



r(r) 



EMti + r)-x) (y(U)-V) 



(5) 



y/T,Mu)-x)* T,Mu)-v) 2 
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Figure 3: Distribution of the average cross-correlation R(t) shown in Fig. |^ and Gaussian 
fit. The standard deviation calculated directly from the histogram coincides with that 
obtained by a Gaussian fit (0.00189 and 0.00191, respectively) 



Let us suppose that GW generate signals in our two detectors with given signal- 
to-noise ratios, say SNR exp i and SNR naut . In the present case, since both detectors have 
about the same sensitivity (T^j 1 = 8.8 mK and T^ff" = 6.1 mK) and are parallel, we take, 
roughly, SNR expi « SNR naut « SNR. We define 

SNRr = ^ (6) 
or 

where or is the standard deviation of our expectations of the dimensionless Rq. This 
quantity can be estimated from the spread of the many available measurements of R(t) 
for different r (see Fig. [J). The distribution of R(r) for all delay times, shown in Fig. |3|, 
can be modeled by a Gaussian (see solid line in the figure), and the standard deviation 
calculated directly from the histogram coincides with that obtained by a Gaussian fit 
(0.00189 and 0.00191, respectively). 

It is possible to evaluate ctr in an alternative way (see Appendix A) , using the 
number of GRB's iV 7 and the number of independent GW data samples N s , in each data 
stretch used for calculating the cross-correlation. The following formula is obtained: 

°R = ; 1 • (7) 

We estimate the number of independent data points by the number of independent data 
of EXPLORER (the bandwidth of NAUTILUS is much narrower, see Table [I]) as follows 

N s = t m Af expl ^ 7200, (8) 

since the detector data are correlated within a time 1/A/. We finally get ctr = 0.0017, in 
excellent agreement with the value 0.0019 deduced from the distribution of R(t), consid- 
ering the roughness of our estimation of the number of independent data, N s . 
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Figure 4: The dimensionless cross-correlation R(t) averaged over the 47 data stretches, 
for various applied signals. In the figures on the left two signals, 100 mK and 50 mK, 
were applied, under the assumption that the GW arrival times be within ±400 s of the 
corresponding GRB arrival times. The figures on the right show the result for applied 
signals of 50 mK and 20 mK, within a time window of ±4s (that is, GW bursts and 
GRB's nearly simultaneous). 

In Appendix A, the following important relationship between the measured quantity 
SNRr (at r = s) and the possible signals expressed in term of signal to noise ratio SNR 
is obtained: 



In order to check the model for the average cross-correlation at zero delay time, and thus 
Eq. |9|, we performed a test by adding to the data, for each GRB, signals of given amplitude 
at the same time for the two detectors. The results are given in Fig.|], where we show the 
cross-correlation for some applied signals for the two time windows of ±400 s and ±4 s. 
The results here are given in terms of the dimensionless average cross-correlation, but 
in the figure we have also indicated the energies of the input signals. We notice that, as 
predicted by Eq. |[ the SNRr increases by reducing the intervals of the cross-correlation, 
for a given SNR. 

3.2 Relations between SNRr, burst energy and dimensionless amplitude 

Since we require the result in terms of GW amplitudes h, we need to relate cross- 
correlation quantities to energies. 

From Eqs. (|J) and @ we obtain 




(9) 




(10) 



Table 2: Comparison of energies calculated, using Eq. fllUp, and simulated (values of an 
input signal that produces SNRr=2). 



window 

s 


E [mK] 
calculated, Eq. fllOD 


E [mK] 
simulated (input signals) 


±400 


43.8 


45 


±40 


24.6 


27 


±16 


19.6 


19 


±4 


13.8 


14 



It would appear that E decreases with increasing N s , i.e. bandwidth, but it can be demon- 
strated that T e ff decreases linearly with increasing bandwidth. 

In terms of the dimensionless amplitude, using Eq. (|3|) with r gw = 1 ms we get, 



h = 7.97 x 10- 18 VE (11) 
with E expressed in kelvin. We finally get, using Eqs. (|6|) and flTOD 




SNR fl = — -\hr (12) 



We observe that a negative value of SNRr gives a negative value for E 2 . To check the 
validity of Eq. fllOD , we applied several signals for various time windows of correlation 
and compared the energies of the input simulated signals to the values calculated using 
Eq. (|mp. This is done in the following way: we note that Rq, in the absence of applied 
signals, is not null for each time window. We determine the signal we must apply in order 
to increase the value of i?o, obtained before the application of the signals, by two standard 
deviations. We then use Eq . (|10"D with the value SNRr = 2. The result is shown in Table |2|. 
The agreement between the values of the simulated input signals and the values calculated 
using Eq. (p~0|) shows that our model is correct. 

4 Inference of the GW burst amplitude 

Having presented the experimental method and the model for the averaged correlation 
at zero delay time Rq, we can infer the values of GW amplitude h consistent with the 
observation. We note that, using Eqs. and H 

energy Eq is related to the measured 

cross-correlation R by 



Hence, the data are summarized by an observed average squared energy Eq = —1.11 x 
10 -3 K 2 , at —1.4 standard deviation from the expected value in the case of noise alone, as 
calculated with the aid of Eq. ( jl2|) where we put SNRr = —1.4. The standard deviation, 
expressed in terms of squared energy, is obtained from Eq. (0), in the case SNRr = 1, 
which gives a E 2 = 0.79 x 10~ 3 K 2 . 



4.1 Probabilistic result and upper limits 

According to the model discussed above, in the case of GW signals of energy E, we expect 
Eq to be a random number, modeled with a Gaussian probability density function around 
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Figure 5: Probability density function f(h\El) (see Eq. The prior used for this 
calculation is a step function. 



E 2 with a standard deviation g E 2\ 



f{E 2 | E 2 ) oc exp 



{El - E 2 ) 2 
2^ 



(14) 



where E is the unknown quantity we wish to infer from the observed value of Eq, given 



in Eq. [13]. This probability inversion is obtained using Bayes' theorem (see, e.g., [22] for 
a physics oriented introduction): 



f(E 2 | El) (x f(E 2 



E 2 )-f (E 2 



(15) 



where f {E 2 ) is the prior probability density function of observing GW signals of squared 
energy E 2 . In fact, we are eventually interested in inferring the GW's amplitude h, related 
to the energy E by Eq. (|TT]). Therefore we have a similar equation: 



f(h\El)<xf(E 2 \h)-f o (h) 



(16) 



where f{El | h) is obtained by a transformation of f(El | E 2 ). As prior for h we considered 
a uniform distribution, bounded to non negative values of h, obtained from Eq. |Tj| i.e. 
fo{h) is a step function 9(h). This seems to us a reasonable choice and it is stable, as 
long as other priors can be conceived which model the "positive attitude of reasonable 
scientists" (see Ref. [[22], [25j) Q. 

The probability density function of h is plotted in Fig. || The highest beliefs are for 
very small values, while values above 1.5 x 1CT 18 are practically ruled out. From Fig. |we 
obtain an expected value and standard deviation for h of 0.56 x 10~ 18 and 0.35 x 10~ 18 , 
respectively, which fully account for what is perceived as a null result. 



A prior distribution alternative to the uniform can be based on the observation that what often seems 
uniform is not the probability per unit of h, but rather the probability per decade of h, i.e. researchers 
may feel equally uncertain about the orders of magnitudes of h. This prior is known as Jeffreys' prior, 
but, in our case, it produces a divergence for h — > in Eq. a direct consequence of the infinite 
orders of magnitudes which are equally believed. To get a finite result we need to put a cut-off at a 
given value of h. This problem is described in depth, for example, in pq| and in pq|. 
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Figure 6: Upper limits for the amplitudes h. Y-axis gives the probability pi{h) < h. 



In these circumstances, we can provide an upper limit, defined as value h(UL), such 
that there is a given probability for the amplitude of GW's to be below it, i.e. 



h(UL) 



f(h | E { 



PL 



(17) 



with pi the chosen probability level. Results are plotted in Fig|]. For example, we can 
exclude the presence of signals of amplitudes h > 1.2 • 10~ 18 with 95% probability. 

4.2 Prior-independent result and upper sensitivity bound 

Probabilistic results depend necessarily on the choice of prior probability density function 
of h. For example, those firmly convinced that GW burst intensities should be in the 
10~ 22 region would never allow a 5% chance to h above 1.2 • 10 -18 . Therefore, in frontier 
research particular care has to be used, before stating probabilistic results. The Bayesian 
approach, thanks to the factorization between likelihood and prior, offers natural ways to 
a present prior-independent result. The simple idea would be just to provide the likelihood 
for each hypothesis under investigation, in our case f{E 2 | h). More conveniently, it has 
been proposed ( ||23|| -p3|) to publish the likelihood rescaled to the asymptotic limit, where 
experimental sensitivity is lost completely; h — 0, in our case. Indicating with 7Z this 
rescaled likelihood, we have 



K{h) 



f(E 2 | h) 



f(E 2 | h ref = 0) 



In statistics jargon, this function gives the Bayes factor of all h hypotheses with respect 
to h = 0. In intuitive terms, it can be interpreted as a "relative belief updating ratio" or 
a " probability density function shape distortion function" , since from Eq. ( |l6l) we have 



f{h\El) 



f{E% | h) 



fo(h) 



f(h = | El) f{El | h = 0) f (h = 0) 
In the present case we get, numerically 

(B-h 4 ) 21 



(19) 



71(h) = A ■ exp 



C 



for h > 



(20) 



10 



R (h) 




Figure 7: Relative belief updating ratio, as a function of the dimensionless amplitude of 
GW's, plotted in log-linear scale. 

where A = 2.6645 is the rescaling factor, B = -4.4 x (10~ 18 ) 4 and C = 20.4 x (KT 18 ) 4 . 

The result is given in Fig. |7], where the choice of the log scale for h is to remember 
that there are infinite orders of magnitudes where the value could be located (and hence 
the problem discussed in footnote [I]). Interpretation of Fig. |7|, in the light of Eqs. flUf)- 
(|T9|) , is straightforward: up to a fraction of 10~ 18 the experimental evidence does not 
produce any change in our belief, while values much larger than 10~ 18 are completely 
ruled out. The region of transition from TZ from 1 to zero identifies a sensitivity bound for 
the experiment. The exact value of this bound is a matter of convention, and could be, for 
example, at TZ = 0.5, or TZ = 0.05. We have 1.3 x 10~ 18 and 1.5 x 10~ 18 , respectively. Note 
that these bounds have no probabilistic meaning. In any case, the full result should be 
considered to be the TZ function, which, being proportional to the likelihood, can easily be 
used to combine results (for independent datasets the global likelihood is the product of 
the likelihoods, and proportional constants can be included in the normalization factor). 
Note that the result given in terms of scaled likelihood and sensitivity bound cannot 
be misleading. In fact, these results are not probabilistic statements about h and no 
one would imagine they were. On the other hand, "confidence limits', which are not 
probabilistic statements on the quantity of interest, tend to be perceived as such (see e.g. 



1 28] and references therein). 



5 Conclusion 

Using for the first time a cross-correlation method applied to the data of two GW detec- 
tors, EXPLORER and NAUTILUS, new experimental upper limits have been determined 
for the burst intensity causing correlations of GW's with GRB's. Analyzing the data over 
47 GRBs, we exclude the presence of signals of amplitude hew > 1-2 • 10~ 18 , with 95% 
probability, with a time window of ± 400 s. With the time window of ± 5 s, we improve 
the previous GW upper limit to about h = 6.5 10~ 19 . 

The result is also given in terms of scaled likelihood and sensitivity bound, which we 
consider the most complete and unbiased way of providing the experimental information. 

In a previous paper we had given more stringent upper limits, but this was 



under the hypothesis that the GW signals always occur at the same time with respect to 
the GRB arrival time. Here, instead, we only require that the time gap between the GRB 
and the GW burst be within a given time window. Similar comparison can be made with 
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the AURIGA/BATSE result |2j|, where an upper limit " h RMS < 1-5 KT 18 with C.L. 
95% " is estimated under the assumption that GW's arrive at the GRB time within a 
time window of ±5 s. 

Finally, we remark that this method can be applied for any expected delay between 
GRB and GW, with appropriate time shifting of the integration window with respect to 
the GBR arrival time, according to the prediction of the chosen model. 
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Appendix A 

Given two independent detectors, let x(t) and y(t) be the measured quantities, which 
are the (filtered) data in our case. 

We introduce the variables £(t) = x(t) — x and r](t) = y(t)—y where x = E[x] = T e ff 
and y = E[y] = T e jf. We recall that the cross- correlation function is 

r(r) = ■ " (t ' + r)) (21) 

with summation extended up to the number of independent samples N s . We calculate 
£ £ 2 = V 2 = N a T? ff . Thus £ 2 • E V 2 = N sTe f f We easily verify that E[t(t)] = 0. 
In the absence of any correlated signal we also have -E^O)] = 0. 

Let us calculate the variance of r(r). We notice that when squaring the numerator 
of Eq.|2T] and taking the average, the cross-terms vanish if the N s data are independent 
from each other. Then, since also £ and r\ are independent variables, we obtain 

, = EMUHU + t)) 2 n s ■ E\ew\ i_ 

{N s T! ff ? (N s T* ff y N s [ ] 

The previous considerations still apply to the cumulative cross-correlation R(r), 
obtained by averaging A 7 independent r(r). The final variance for the cross-correlation 
R(t) is 

^ = (23) 

Let us now consider an energy signal S on both detectors at the same time. The 
expected value of the cross-correlation r(r) at r = will be positive. In the case of N s 
independent data the signal S will appear in one datum only and we have: 

E[R(r = 0) - R(r + 0)] = TF7FS - = (24) 



where we put SNR = ^Ay. We note that, in Eq.|24], E[R(t ^ 0)] = also when a signal 
is present. 

Putting also 

SN R R = EMl^M (25 ) 
(as already been defined in the text in Eq. |6]) we obtain 

SNR fl = SNR 2 A /5- (26) 
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